(*
  Copyright (c) 2009 Barry Schwartz

  Permission is hereby granted, free of charge, to any person
  obtaining a copy of this software and associated documentation
  files (the "Software"), to deal in the Software without
  restriction, including without limitation the rights to use,
  copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the
  Software is furnished to do so, subject to the following
  conditions:

  The above copyright notice and this permission notice shall be
  included in all copies or substantial portions of the Software.

  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  OTHER DEALINGS IN THE SOFTWARE.
*)

open XNum
  
type spline = (num * num * num * num * num * num * num * num)

type side_spec =
  | Open
  | Curl of num
  | Angle of float
  | Explicit of num * num
  | Endpoint
  | Cycle

type knot =
  { left_spec : side_spec; right_spec : side_spec; left_tension : num;
    right_tension : num; point_x : num; point_y : num
  }

type path_spec = knot list

let log_knot k =
  (Printf.eprintf "(%f,%f)" (float_of_num k.point_x) (float_of_num k.point_y);
   (match k.left_spec with
    | Open -> Printf.eprintf " open"
    | Curl c -> Printf.eprintf " curl %f" (float_of_num c)
    | Angle a -> Printf.eprintf " angle %f" a
    | Explicit (x, y) ->
        Printf.eprintf " explicit %f %f" (float_of_num x) (float_of_num y)
    | Endpoint -> Printf.eprintf " endpoint"
    | Cycle -> Printf.eprintf " cycle");
   Printf.eprintf " %f" (float_of_num k.left_tension);
   (match k.right_spec with
    | Open -> Printf.eprintf " open"
    | Curl c -> Printf.eprintf " curl %f" (float_of_num c)
    | Angle a -> Printf.eprintf " angle %f" a
    | Explicit (x, y) ->
        Printf.eprintf " explicit %f %f" (float_of_num x) (float_of_num y)
    | Endpoint -> Printf.eprintf " endpoint"
    | Cycle -> Printf.eprintf " cycle");
   Printf.eprintf " %f\n" (float_of_num k.right_tension))
  
let log_knots knots =
  for i = 0 to (Array.length knots) - 1 do
    Printf.eprintf "%d: " i;
    log_knot knots.(i)
  done
  
(* side specifications *)
let is_open spec = match spec with | Open -> true | _ -> false
  
let is_curl spec = match spec with | Curl _ -> true | _ -> false
  
let is_angle spec = match spec with | Angle _ -> true | _ -> false
  
let is_explicit spec = match spec with | Explicit (_, _) -> true | _ -> false
  
let is_endpoint spec = match spec with | Endpoint -> true | _ -> false
  
let is_cycle spec = match spec with | Cycle -> true | _ -> false
  
let curl spec =
  match spec with | Curl c -> c | _ -> invalid_arg "Bezier.curl"
  
let convert_open_to_curl spec =
  match spec with | Open -> Curl num_one | _ -> spec
  
(* geometric functions *)
let pi = 4.0 *. (atan 1.0)
  
let angle_of_vec x y =
  (180.0 /. pi) *. (atan2 (float_of_num y) (float_of_num x))
  
let sind x = sin ((pi *. x) /. 180.0)
  
let cosd x = cos ((pi *. x) /. 180.0)
  
let fast_reduce_angle a =
  if a > 180.0 then a -. 180.0 else if a < (-180.0) then a +. 180.0 else a
  
let slow_reduce_angle a = mod_float a 360.0
  
let num_4 = num_of_int 4
  
let num_16 = num_of_int 16
  
let num_sqrt_2 = num_of_float (sqrt 2.0)
  
let num_sqrt_5 = num_of_float (sqrt 5.0)
  
let const1 = mult_num (num_of_ints 3 2) (sub_num num_sqrt_5 num_one)
  
let const2 = mult_num (num_of_ints 3 2) (sub_num num_three num_sqrt_5)
  
let velocity sin_theta cos_theta sin_phi cos_phi tension =
  let n =
    add_num num_two
      (mult_num
         (mult_num
            (mult_num num_sqrt_2 (sub_num sin_theta (div_num sin_phi num_16)))
            (sub_num sin_phi (div_num sin_theta num_16)))
         (sub_num cos_theta cos_phi)) in
  let d =
    add_num (add_num num_three (mult_num const1 cos_theta))
      (mult_num const2 cos_phi) in
  let x = div_num n (mult_num d tension)
  in if ge_num x num_4 then num_4 else x
  
let curl_ratio curl left_tension right_tension =
  let alpha = div_num num_one right_tension in
  let beta = div_num num_one left_tension in
  let r = div_num alpha beta in
  let gamma = mult_num (mult_num r r) curl in
  let d = sub_num (add_num (mult_num gamma alpha) num_three) beta in
  let n = add_num (mult_num gamma (sub_num num_three alpha)) beta in
  let x = div_num n d in if ge_num x num_4 then num_4 else x
  
(* Join knots with the same coordinates by a straight line. *)
let join_points spec =
  let len = Array.length spec
  in
    for i = 0 to len - 1 do
      let n = if (i + 1) < len then i + 1 else 0 in
      let k1 = spec.(i) in
      let k2 = spec.(n)
      in
        if (eq_num k1.point_x k2.point_x) && (eq_num k1.point_y k2.point_y)
        then
          (* add straight line (of length 0!) between <k> and <last_knot> *)
          (match k1.right_spec with
           | Open | Curl _ | Angle _ ->
               (spec.(i) <-
                  {
                    (k1)
                    with
                    right_spec = Explicit (k2.point_x, k2.point_y);
                    left_spec = convert_open_to_curl k1.left_spec;
                  };
                spec.(n) <-
                  {
                    (k2)
                    with
                    left_spec = Explicit (k1.point_x, k1.point_y);
                    right_spec = convert_open_to_curl k2.right_spec;
                  })
           | _ -> ())
        else ()
    done
  
let calculate_angles spec =
  let len = Array.length spec in
  let delta_x = Array.make len num_zero in
  (* vector from one point to the next one *)
  let delta_y = Array.make len num_zero in
  let dist = Array.make len num_zero in
  (* distance between the current point and the next one *)
  let psi = Array.make len num_zero
  in
    ((* turning angle at the corresponding point *)
     for i = 0 to len - 1 do
       (let n = if (i + 1) < len then i + 1 else 0 in
        let k1 = spec.(i) in
        let k2 = spec.(n)
        in
          (delta_x.(i) <- sub_num k2.point_x k1.point_x;
           delta_y.(i) <- sub_num k2.point_y k1.point_y;
           dist.(i) <-
             num_of_float
               (sqrt
                  (float_of_num
                     (add_num (mult_num delta_x.(i) delta_x.(i))
                        (mult_num delta_y.(i) delta_y.(i)))));
           if i > 0
           then
             (let sine = div_num delta_y.(i - 1) dist.(i - 1) in
              let cosine = div_num delta_x.(i - 1) dist.(i - 1)
              in
                psi.(i) <-
                  num_of_float
                    (angle_of_vec
                       (add_num (mult_num delta_x.(i) cosine)
                          (mult_num delta_y.(i) sine))
                       (sub_num (mult_num delta_y.(i) cosine)
                          (mult_num delta_x.(i) sine))))
           else ()))
     done;
     if is_cycle spec.(0).left_spec
     then
       (let sine = div_num delta_y.(len - 1) dist.(len - 1) in
        let cosine = div_num delta_x.(len - 1) dist.(len - 1)
        in
          (psi.(0) <-
             num_of_float
               (angle_of_vec
                  (add_num (mult_num delta_x.(0) cosine)
                     (mult_num delta_y.(0) sine))
                  (sub_num (mult_num delta_y.(0) cosine)
                     (mult_num delta_x.(0) sine)));
           psi.(len - 1) <- psi.(0)))
     else (psi.(0) <- num_zero; psi.(len - 1) <- num_zero);
     (delta_x, delta_y, dist, psi))
  
let remove_open_at_endpoints spec first_knot last_knot =
  let k1 = spec.(first_knot) in
  let k2 = spec.(last_knot)
  in
    (if is_open k1.right_spec
     then
       (match k1.left_spec with
        | Explicit (x, y) ->
            let dx = sub_num k1.point_x x in
            let dy = sub_num k1.point_y y
            in
              if (eq_num dx num_zero) && (eq_num dy num_zero)
              then
                spec.(first_knot) <- { (k1) with right_spec = Curl num_one; }
              else
                spec.(first_knot) <-
                  { (k1) with right_spec = Angle (angle_of_vec dx dy); }
        | _ -> ())
     else ();
     if is_open k2.left_spec
     then
       (match k2.right_spec with
        | Explicit (x, y) ->
            let dx = sub_num x k2.point_x in
            let dy = sub_num y k2.point_y
            in
              if (eq_num dx num_zero) && (eq_num dy num_zero)
              then
                spec.(last_knot) <- { (k2) with left_spec = Curl num_one; }
              else
                spec.(last_knot) <-
                  { (k2) with left_spec = Angle (angle_of_vec dx dy); }
        | _ -> ())
     else ())
  
let set_control_points knots left_knot right_knot delta_x delta_y sin_theta
                       cos_theta sin_phi cos_phi =
  let lk = knots.(left_knot) in
  let rk = knots.(right_knot) in
  let left_tension = lk.right_tension in
  let right_tension = rk.left_tension in
  let left_vel =
    velocity sin_theta cos_theta sin_phi cos_phi (abs_num left_tension) in
  let right_vel =
    velocity sin_phi cos_phi sin_theta cos_theta (abs_num right_tension) in
  let (left_vel2, right_vel2) =
    if
      ((lt_num left_tension num_zero) || (lt_num right_tension num_zero)) &&
        (((ge_num sin_theta num_zero) && (ge_num sin_phi num_zero)) ||
           ((le_num sin_theta num_zero) && (le_num sin_phi num_zero)))
    then
      (let st = abs_num sin_theta in
       let sp = abs_num sin_phi in
       let s =
         mult_num (num_of_ints 4097 4096)
           (add_num (mult_num st cos_phi) (mult_num sp cos_theta))
       in
         if gt_num s num_zero
         then
           (let l =
              if lt_num left_tension num_zero
              then min_num (div_num sp s) left_vel
              else left_vel in
            let r =
              if lt_num right_tension num_zero
              then min_num (div_num st s) right_vel
              else right_vel
            in (l, r))
         else (left_vel, right_vel))
    else (left_vel, right_vel) in
  let left_x =
    add_num lk.point_x
      (mult_num left_vel2
         (sub_num (mult_num delta_x cos_theta) (mult_num delta_y sin_theta))) in
  let left_y =
    add_num lk.point_y
      (mult_num left_vel2
         (add_num (mult_num delta_y cos_theta) (mult_num delta_x sin_theta))) in
  let right_x =
    sub_num rk.point_x
      (mult_num right_vel2
         (add_num (mult_num delta_x cos_phi) (mult_num delta_y sin_phi))) in
  let right_y =
    sub_num rk.point_y
      (mult_num right_vel2
         (sub_num (mult_num delta_y cos_phi) (mult_num delta_x sin_phi)))
  in
    (knots.(left_knot) <-
       { (lk) with right_spec = Explicit (left_x, left_y); };
     knots.(right_knot) <-
       { (rk) with left_spec = Explicit (right_x, right_y); })
  
let calculate_splines knots (delta_x, delta_y, dist, psi) first_knot
                      last_knot =
  (* We solve a system of equations of the form  x_i + u_i x_{i+1} = v_i + w_i x_0 *)
  let len = Array.length knots in
  let idx i = if i < 0 then i + len else if i >= len then i - len else i in
  let u = Array.make len num_zero in
  let v = Array.make len num_zero in
  let w = Array.make len num_zero in
  let theta = Array.make len num_zero in
  let next = idx (first_knot + 1)
  in
    ((match knots.(first_knot).right_spec with
      | Angle a ->
          (match knots.(next).left_spec with
           | Angle a2 ->
               let angle =
                 angle_of_vec delta_x.(first_knot) delta_y.(first_knot)
               in
                 set_control_points knots first_knot next
                   delta_x.(first_knot) delta_y.(first_knot)
                   (num_of_float (sind (a -. angle)))
                   (num_of_float (cosd (a -. angle)))
                   (num_of_float (sind (a2 -. angle)))
                   (num_of_float (cosd (a2 -. angle)))
           | _ ->
               ((* theta_0 + 0 * theta_1 = a + 0 * theta_0 *)
                u.(first_knot) <- num_zero;
                v.(first_knot) <-
                  num_of_float
                    (fast_reduce_angle
                       (a -.
                          (angle_of_vec delta_x.(first_knot)
                             delta_y.(first_knot))));
                w.(first_knot) <- num_zero))
      | Curl c ->
          (match knots.(next).left_spec with
           | Curl _ ->
               let kl = knots.(first_knot) in
               let kr = knots.(next) in
               let sl = mult_num num_three (abs_num kl.right_tension) in
               let sr = mult_num num_three (abs_num kr.left_tension)
               in
                 (knots.(first_knot) <-
                    {
                      (kl)
                      with
                      right_spec =
                        Explicit
                          (add_num kl.point_x
                             (div_num delta_x.(first_knot) sl),
                          add_num kl.point_y
                            (div_num delta_y.(first_knot) sl));
                    };
                  knots.(next) <-
                    {
                      (kr)
                      with
                      left_spec =
                        Explicit
                          (sub_num kr.point_x
                             (div_num delta_x.(first_knot) sr),
                          sub_num kr.point_y
                            (div_num delta_y.(first_knot) sr));
                    })
           | _ ->
               let lt = abs_num knots.(first_knot).right_tension in
               let rt = abs_num knots.(next).left_tension
               in
                 (if (eq_num lt num_one) && (eq_num rt num_one)
                  then
                    u.(first_knot) <-
                      div_num (add_num (add_num c c) num_one)
                        (add_num c num_two)
                  else u.(first_knot) <- curl_ratio c lt rt;
                  v.(first_knot) <-
                    minus_num (mult_num psi.(next) u.(first_knot));
                  w.(first_knot) <- num_zero))
      | Open ->
          (u.(first_knot) <- num_zero;
           v.(first_knot) <- num_zero;
           w.(first_knot) <- num_one)
      | _ -> assert false);
     (let rec iter middle =
        let left = idx (middle - 1) in
        let right = idx (middle + 1)
        in
          ((match knots.(middle).left_spec with
            | Cycle | Open ->
                let lt = abs_num knots.(left).right_tension in
                let rt = abs_num knots.(right).left_tension in
                let a =
                  div_num num_one (sub_num (mult_num num_three lt) num_one) in
                let b =
                  div_num num_one (sub_num (mult_num num_three rt) num_one) in
                let c = sub_num num_one (mult_num a u.(left)) in
                let d =
                  mult_num dist.(middle)
                    (sub_num num_three (div_num num_one lt)) in
                let e =
                  mult_num dist.(left)
                    (sub_num num_three (div_num num_one rt)) in
                let s =
                  abs_num
                    (div_num knots.(middle).left_tension
                       knots.(middle).right_tension) in
                let f =
                  div_num e
                    (add_num e (mult_num (mult_num (mult_num c d) s) s))
                in
                  (u.(middle) <- mult_num b f;
                   let acc = minus_num (mult_num psi.(right) u.(middle))
                   in
                     (if is_curl knots.(left).right_spec
                      then
                        (v.(middle) <-
                           sub_num acc
                             (mult_num psi.(idx (first_knot + 1))
                                (sub_num num_one f));
                         w.(middle) <- num_zero)
                      else
                        (let f = div_num (sub_num num_one f) c in
                         let acc = sub_num acc (mult_num psi.(middle) f) in
                         let f = mult_num a f
                         in
                           (v.(middle) <- sub_num acc (mult_num f v.(left));
                            w.(middle) <- minus_num (mult_num f w.(left))));
                      if is_cycle knots.(middle).left_spec
                      then (* set theta.(first_knot) = theta.(last_knot) *)
                        (let rec iter2 a b i =
                           let c = sub_num v.(i) (mult_num a u.(i)) in
                           let d = sub_num w.(i) (mult_num b u.(i))
                           in
                             if i <> last_knot
                             then iter2 c d (idx (i - 1))
                             else
                               (let e = div_num c (sub_num num_one d)
                                in
                                  (theta.(last_knot) <- e;
                                   v.(first_knot) <- e;
                                   let rec iter3 i =
                                     let n = idx (i + 1)
                                     in
                                       (v.(i) <-
                                          add_num v.(i) (mult_num e w.(i));
                                        if n <> last_knot
                                        then iter3 n
                                        else ())
                                   in iter3 (idx (first_knot + 1))))
                         in iter2 num_zero num_one left)
                      else ()))
            | Curl c ->
                let lt = abs_num knots.(middle).left_tension in
                let rt = abs_num knots.(left).right_tension in
                let s =
                  if (eq_num lt num_one) && (eq_num rt num_one)
                  then
                    div_num (add_num (add_num c c) num_one)
                      (add_num c num_two)
                  else curl_ratio c rt lt
                in
                  theta.(middle) <-
                    div_num (minus_num (mult_num s v.(left)))
                      (sub_num num_one (mult_num s u.(left)))
            | Angle a ->
                theta.(middle) <-
                  num_of_float
                    (fast_reduce_angle
                       (a -. (angle_of_vec delta_x.(left) delta_y.(left))))
            | _ -> ());
           (* We have already computed the values. *)
           if middle <> last_knot then iter right else ())
      in iter (idx (first_knot + 1)));
     let rec iter i =
       let n = idx (i + 1)
       in
         (if i <> last_knot
          then (* If the path is cyclic then we have already        *)
            theta.(i) <- sub_num v.(i) (mult_num u.(i) theta.(n))
          else (* computed theta.(i) for i = first_knot = last_knot *) ();
          let phi = sub_num (minus_num psi.(n)) theta.(n) in
          let sin_theta = sind (float_of_num theta.(i)) in
          let cos_theta = cosd (float_of_num theta.(i)) in
          let sin_phi = sind (float_of_num phi) in
          let cos_phi = cosd (float_of_num phi)
          in
            (set_control_points knots i n delta_x.(i) delta_y.(i)
               (num_of_float sin_theta) (num_of_float cos_theta)
               (num_of_float sin_phi) (num_of_float cos_phi);
             if i <> first_knot then iter (idx (i - 1)) else ()))
     in iter (idx (last_knot - 1)))
  
let compute_path spec =
  let knots = Array.of_list spec in
  let len = Array.length knots
  in
    if len < 2
    then
      if len = 0
      then [|  |]
      else
        (let k = knots.(0)
         in
           [|
             ((k.point_x), (k.point_y), (k.point_x), (k.point_y),
              (k.point_x), (k.point_y), (k.point_x), (k.point_y))
           |])
    else
      (join_points knots;
       (let dist_psi = calculate_angles knots in
        let idx i =
          if i < 0 then i + len else if i >= len then i - len else i in
        let make_path knots =
          let n =
            if is_endpoint knots.(0).left_spec
            then (Array.length knots) - 1
            else Array.length knots in
          let make_spline i =
            let k1 = knots.(i) in
            let k2 = knots.(idx (i + 1))
            in
              match ((k1.right_spec), (k2.left_spec)) with
              | (Explicit (x1, y1), Explicit (x2, y2)) ->
                  ((k1.point_x), (k1.point_y), x1, y1, x2, y2, (k2.point_x),
                   (k2.point_y))
              | _ -> assert false
          in Array.init n make_spline in
        let rec find_next_break i =
          match ((knots.(i).left_spec), (knots.(i).right_spec)) with
          | (Open, Open) -> find_next_break (idx (i + 1))
          | _ -> i in
        let first_break =
          match knots.(0).left_spec with
          | Endpoint -> 0
          | Cycle -> find_next_break 1
          | _ -> assert false
        in
          (if first_break > 0
           then
             (* The path is cyclic and it contains a real break point.
         We replace the |Cycle| by |Open|. *)
             knots.(0) <- { (knots.(0)) with left_spec = Open; }
           else ();
           let rec iter cur_break =
             let continue next_break =
               if
                 (next_break <> first_break) &&
                   (not (is_endpoint knots.(next_break).right_spec))
               then iter next_break
               else make_path knots
             in
               if is_explicit knots.(cur_break).right_spec
               then (* nothing to do *) continue (idx (cur_break + 1))
               else
                 (let next_break = find_next_break (idx (cur_break + 1))
                  in
                    (remove_open_at_endpoints knots cur_break next_break;
                     calculate_splines knots dist_psi cur_break next_break;
                     continue next_break))
           in iter first_break)))
  
(* building path specifications *)
let make_spec x y =
  [ {
      left_spec = Endpoint;
      right_spec = Open;
      left_tension = num_one;
      right_tension = num_one;
      point_x = x;
      point_y = y;
    } ]
  
let make_specs_consistent knot =
  match ((knot.left_spec), (knot.right_spec)) with
  | (Open, Curl c) -> { (knot) with left_spec = Curl c; }
  | (Open, Angle a) -> { (knot) with left_spec = Angle a; }
  | (Open, Endpoint) -> { (knot) with left_spec = Curl num_one; }
  | (Curl c, Open) -> { (knot) with right_spec = Curl c; }
  | (Angle a, Open) -> { (knot) with right_spec = Angle a; }
  | (Endpoint, Open) -> { (knot) with right_spec = Curl num_one; }
  | _ -> knot
  
let rec cleanup_spec cycle result spec =
  match spec with
  | [ k ] ->
      if cycle
      then
        (let x = make_specs_consistent { (k) with left_spec = Open; } in
         let y = { (x) with left_spec = Cycle; } in y :: result)
      else (make_specs_consistent k) :: result
  | k :: ks -> cleanup_spec cycle ((make_specs_consistent k) :: result) ks
  | [] -> assert false
  
let close_spec spec cycle =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_cycle k.right_spec
      then invalid_arg "Bezier.close_spec: last point misssing!"
      else
        (let s =
           if cycle then spec else { (k) with right_spec = Endpoint; } :: ks
         in compute_path (cleanup_spec cycle [] s))
  
let add_point spec x y =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_cycle k.right_spec
      then (* |Cycle| indicates that the coordinates of <k> are invalid *)
        { (k) with right_spec = Open; point_x = x; point_y = y; } :: ks
      else
        {
          left_spec = Open;
          right_spec = Open;
          left_tension = num_one;
          right_tension = num_one;
          point_x = x;
          point_y = y;
        } :: spec
  
let add_in_dir spec angle =
  match spec with
  | [] -> assert false
  | k :: _ ->
      if is_cycle k.right_spec
      then invalid_arg "Bezier.add_in_dir: left side already specified!"
      else
        {
          left_spec = Angle (slow_reduce_angle angle);
          right_spec = Cycle;
          left_tension = num_one;
          right_tension = num_one;
          point_x = num_zero;
          point_y = num_zero;
        } :: spec
  
let add_in_curl spec curl =
  match spec with
  | [] -> assert false
  | k :: _ ->
      if is_cycle k.right_spec
      then invalid_arg "Bezier.add_in_curl: left side already specified!"
      else
        {
          left_spec = Curl curl;
          right_spec = Cycle;
          left_tension = num_one;
          right_tension = num_one;
          point_x = num_zero;
          point_y = num_zero;
        } :: spec
  
let add_in_tension spec tension =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_cycle k.right_spec
      then { (k) with left_tension = tension; } :: ks
      else
        {
          left_spec = Open;
          right_spec = Cycle;
          left_tension = tension;
          right_tension = num_one;
          point_x = num_zero;
          point_y = num_zero;
        } :: spec
  
let add_out_dir spec angle =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_open k.right_spec
      then { (k) with right_spec = Angle (slow_reduce_angle angle); } :: ks
      else invalid_arg "Bezier.add_out_dir: right side already specified!"
  
let add_out_curl spec curl =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_open k.right_spec
      then { (k) with right_spec = Curl curl; } :: ks
      else invalid_arg "Bezier.add_out_curl: right side already specified!"
  
let add_out_tension spec tension =
  match spec with
  | [] -> assert false
  | k :: ks -> { (k) with right_tension = tension; } :: ks
  
let add_control_points spec x1 y1 x2 y2 =
  match spec with
  | [] -> assert false
  | k :: ks ->
      if is_open k.right_spec
      then
        {
          left_spec = Explicit (x2, y2);
          right_spec = Cycle;
          left_tension = num_one;
          right_tension = num_one;
          point_x = num_zero;
          point_y = num_zero;
        } :: { (k) with right_spec = Explicit (x1, y1); } :: ks
      else
        if is_cycle k.right_spec
        then invalid_arg "Bezier.add_control_points: point missing!"
        else
          invalid_arg
            "Bezier.add_control_points: right side already specified!"
  

